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ABSTRACT 


The  Infrasound  Experts  Group  of  the  Geneva  Conference  on  Disarmament  Ad  Hoc  Committee  on 
a  Nuclear  Test  Ban  has  recommended  an  infrasound  array  design  consisting  of  four  elements, 
with  three  elements  forming  an  equilateral  triangle  and  the  fourth  at  the  center  of  the  triangle. 

The  Experts  recommended  that  the  sides  of  the  triangle  be  in  the  range  of  1  to  3  kilometers  (km). 

In  this  report,  correlation  as  a  function  of  period  and  sensor  spacing,  and  signal-to-noise  ratios 
(S/N)  as  a  function  of  period  are  estimated  from  atmospheric  data  from  nuclear  explosions  of  2.2 
and  2.6  kilotons  (kt)  recorded  at  sensors  with  spacings  near  1  km  and  at  distances  from  the  explo¬ 
sions  of  700  km.  The  correlation  estimates  are  found  to  be  consistent  with  parameters  describing 
the  spread  of  the  signal  in  wavenumber  space  as  discussed  by  Blandford  (1997).  These  parame¬ 
ters  were  first  estimated  by  Mack  and  Flinn  (1971)  for  events  of  much  larger  yield  at  larger  dis¬ 
tances,  longer  period,  and  greater  inter-sensor  spacings.  Data  at  10  km  spacing  from  a  2.2  kt 
event  at  a  distance  of 2200  km  is  found  to  be  consistent  with  the  Mack  and  Flinn  parameters.  The 
correlation  and  S/N  estimates  are  used  to  show  that,  for  most  International  Monitoring  Commu¬ 
nity  applications,  the  optimum  period  for  detection  and  azimuth  estimation  is  1  second,  and  that 
the  optimum  array  aperture  is  1  km. 

Preliminary  analysis  of  a  signal  at  the  TXAR  infrasonic  array  from  a  White  Sands  10-ton  chemi¬ 
cal  explosion  at  a  distance  of 450  km  gives  much  higher  signal  correlations  at  1  Hertz  (Hz)  and 
1  km  spacing,  suggesting  that  one-hop  signals  have  approximately  1/5  the  dispersion  in  azimuth 
and  velocity  as  do  multi-hop  signals.  Thus,  paradoxically,  it  may  be  that  an  array  larger  than  1  km 
may  be  optimum  for  events  at  distances  <500  km,  while  being  suboptimal  for  events  at  greater 
distances. 
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INTRODUCTION 


In  an  earlier  report  on  infrasound  array  design,  Blandford  (1997)  used  coherence  functions  deter¬ 
mined  by  Mack  and  Flinn  (1971)  (hereafter  referred  to  as  MF)  from  infrasound  wave  data  from 
one  of  China’s  megaton-range  explosions  recorded  at  the  Montana  Large  Aperture  Meteorologi¬ 
cal  Array  (LAMA).  The  coherence  was  determined  using  periods  T  >10  seconds  (sec),  and  at 
sensor  intervals  of  x  >7  kilometers  (km).  Using  these  results,  Blandford  (1997)  showed  that  for  a 
“detection”  beam  root-mean-square  (rms)  signal-to-noise  (S/N)  amplitude  ratio  of  1.5  at  the 
expected  detection  period  of  5  seconds  for  a  1  kiloton  (kt)  explosion,  a  1-km  aperture  array  may 
be  expected  to  have  an  azimuth  estimation  standard  error  smaller  than  the  experimentally 
observed  1.8°  true  azimuthal  error.  Thus,  it  appeared  that  an  array  aperture  larger  than  1  km, 
which  would  likely  entail  increased  logistic  expenses,  would  not  be  necessary  for  adequate  azi¬ 
muth  estimation  capability. 

In  addition,  it  was  shown,  assuming  the  MF  data  parameterization  could  be  used  to  extrapolate  to 
1-sec  period,  that  for  array  apertures  of  2-3  km,  1  Hertz  (Hz)  signals  would  be  substantially 
uncorrelated  so  that  array  signal  estimates  at  this  period  would  be  poor.  This  could  be  undesirable 
if  1  Hz  data  were  to  be  used  for  discrimination  between  explosion  signals  and  other  signals.  Thus, 
again,  a  1-km  aperture  was  indicated. 

A  weakness  in  the  work  of  Blandford  (1997)  was  that  there  were  no  data  to  show  that  the  MF  data 
could  be  reliably  extrapolated  to  1  Hz.  Moreover,  there  were  no  data  available  on  the  variation  of 
S/N  as  a  function  of  period  from  sources  of  interest.  There  were,  however,  data  cited  to  show  that 
the  peak  in  the  signal  spectrum  from  1  kt  events  at  1000  km  is  near  5  sec  period,  but  it  is  not  nec¬ 
essarily  the  case  that  this  is  the  same  as  the  peak  in  S/N,  which  is  the  critical  parameter  for  detec¬ 
tion  and  azimuth  estimation.  In  addition,  for  both  detection  and  azimuth  estimation,  capability 
increases  in  proportion  to  the  time-bandwidth  product.  Thus,  optimum  capability  may  occur  at  a 
frequency  higher  than  that  for  which  there  is  the  highest  S/N. 

In  this  report,  0.5-  to  10-sec  signals  from  the  2.6  kt  nuclear  atmospheric  test,  Tanana,  recorded  by 
a  1-km  aperture  array  on  Palmyra  Island  are  examined.  Based  on  these  and  other  similar  data,  we 
find  that  extrapolation  of  the  MF  functions  is  justified.  Consequently,  the  interaction  of  observed 
S/N  with  array  aperture  implies  that  the  best  1-kt  explosion  detection  and  azimuth  estimation  per¬ 
formance  is  at  1-sec  period  and  at  an  aperture  of  1  km. 

Figure  1,  from  Blandford  (1997),  illustrates  the  squared  coherence  as  a  function  of  period  and  ele¬ 
ment  spacing,  both  parallel  and  perpendicular  to  the  wavefront  calculated  from  equation  (1). 
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In  equation  (1)  y  is  coherence,  k0  is  the  wavenumber,  A0  is  the  spread  in  azimuth  of  the  signal 
packet,  M  is  the  spread  in  wavenumber,  and  x  and  y  are  the  relative  locations  of  the  two  sensors. 
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Infrasound  Squared  Coherence,  Parallel . and  Perpendicular 


Figure  1.  Squared  coherence  as  a  function  of  intersensor  distance  for  T=10,  5,  2,  and  1 
seconds  according  to  the  equation  above,  using  Ac=0.015  km/sec  and  A©=5°  as 
determined  from  Mack  and  Flinn  (1971),  parallel  and  perpendicular  to  the  wavefront. 
Observed  correlation  estimates  should  be  compared  to  the  theoretical  estimates  of 
coherence  and  not  to  squared  coherence. 


In  subsequent  sections,  we  discuss  the  Palmyra  data,  present  the  methods  used  to  determine  the 
Palmyra  signal  correlation  values  as  a  function  of  period,  correcting  for  the  noise  (which  is  a 
problem  at  low  S/N),  compare  the  correlation  results  to  the  MF  results,  develop  an  “F  detector” 
detection  theory  with  which  to  compare  the  detection  capability  of  different  array  apertures,  and 
determine  the  optimum  array  aperture  with  respect  to  both  detection  and  azimuth  estimation 
(using  the  azimuth  estimation  theory  discussed  in  Blandford  (1997)).  This  theory,  by  Shumway  et 
al.  (1999),  is  also  given  in  this  report  in  Appendix  1,  and  is  extended  there  to  include  results  for 
simple  beamforming  as  well  as  for  the  case  of  optimum  beamforming,  which  was  given  in 
Blandford  (1997).  We  then  show  that  the  same  aperture  results  are  obtained  if  simple  beamform¬ 
ing  is  used  for  azimuth  estimation  in  place  of  optimum  (weighted)  beamforming,  and  also  present 
some  additional  coherence  data  which  also  is  consistent  with  MF.  The  paper  closes  with  a  sum¬ 
mary  and  conclusions. 
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PALMYRA  AND  OTHER  SMALL  INFRASOUND  ARRAYS, 

AND  PALMYRA  DATA 


The  2.6  kt  atmospheric  event  Tanana  was  detonated  at  1.65°N,  157.28°W  at  16:08:52  on  25  May 
1962  as  part  of  Operation  Dominic.  The  shot  point  was  south  of  Kiritimati  Island,  1.87°N, 
157.33°W  (at  shot  time  the  island  was  a  United  Kingdom  possession  and  was  named  Christmas 
Island). 

Approximately  700  km  to  the  northwest,  on  Palmyra  Island,  a  US  possession  at  5.88°N, 
162.08°W,  the  US  Government  in  1962  installed  die  temporary  infrasound  station  “Small  Fry.” 
As  seen  in  Figure  2,  Palmyra  is  made  up  of  a  number  of  islands  with  a  typical  dimension  of  1  km, 
connected  by  several  man-made  structures,  erected  just  before  WWII,  within  a  reef  approximately 
3  km  by  10  km. 

Thus,  although  the  original  coordinates  of  the  1962, 3 -element  Small  Fry  array  have  not  yet  been 
recovered  from  archives,  it  seems  plausible  that  the  elements  could  be  spaced  at  distances  on  the 
order  of  1  km.  As  we  shall  see,  this  is  supported  by  the  observation  that  the  maximum  signal 
delay  across  the  array  is  approximately  3  seconds.  The  existence  of  this  small  array  is  critical 
because  data  from  it  comprises  almost  all  the  atmospheric  nuclear  explosion  data  that  is  available 
at  the  lower  limit  of  instrument  spacings  which  are  representative  of  the  proposed  infrasound 
arrays  to  be  built  for  International  Monitoring  Community  monitoring. 

Perhaps  unfortunately  for  our  present  purposes,  the  Small  Fry  array  may  have  had  linear  pipe  fil¬ 
ters  approximately  300  meters  in  length;  it  is  believed  that  such  filters  were  designed  and  became 
standard  in  US  Government  arrays  soon  after  1955  since  the  principal  objective  of  the  network 
was  to  detect  1 00  kt  or  greater  explosions  at  distances  >3000  kilometers.  This  length  of  pipe  array 
would  be  about  equal  to  the  wavelength  of  a  1  Hz  signal.  The  signal  at  1  Hz  would,  therefore, 
likely  be  strongly  attenuated  along  some  azimuths,  thus  leading  to  an  underestimate  of  S/N  at/>l 
Hz  and,  therefore,  to  an  underestimation  of  the  capability  of  1  Hz  detection  and  azimuth  estima¬ 
tion  for  arrays  of  instruments  with  pipe  arrays  of,  say,  100  meters  characteristic  dimension. 

An  initial  infrasound  array  was  installed  on  Palmyra  in  1951  for  Operation  Greenhouse.  A  Febru¬ 
ary  1952  Final  Report  of  the  “Ascension  Group”  is  the  source  of  Figure  2.  Reference  to  a  naviga¬ 
tion  chart  for  Palmyra  suggests  that  the  four  elements  of  the  1952  array  were  located  at  sites 
which  can  be  described  as: 

(1)  Strawn  Island 

(2)  Aviation  Island,  0.25  km  East  of  the  North  end  of  the  North-South  Causeway 

(3)  NE  end  of  Eastern  Island 

(4)  Eastern  end  of  Kaula  Island,  0.2  km  West  of  the  South  end  of  the  North-South  Causeway. 

For  this  earlier  array,  the  greatest  separation  was  between  elements  1  and  4,  and  was  4.5  km;  and 
the  least  was  between  elements  2  and  3,  and  was  2  km.  Thus,  it  is  apparent  that  the  1962  Small 
Fry  array  must  have  been  deployed  differently  from  the  1951  array.  There  were  also  other  small- 
to-intermediate  size  arrays  deployed  in  1951  at  the  Majuro  Atoll,  and  at  Eniwetok  itself  near 
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which  the  four  explosions  of  operation  Greenhouse  were  detonated.  The  1952  data  are  currently 
being  recovered  and  digitized. 

Figure  3  shows  the  unfiltered  waveforms  of  Tanana  at  Palmyra  with  “A-B”  and  “C-D”  indicating 
two  100-second  time  windows  in  the  signals  which  are  used  for  analysis,  and  Figure  4  shows 
channel  3  filtered  through  a  set  of  passbands  centered  at  the  periods  indicated. 

The  Tanana  signal  in  Figure  3  is  somewhat  misleading;  it  might  seem  that  the  time  axis  has  been 
reversed,  that  the  earliest  times  are  on  the  right  of  the  plot  and  that  the  signal  first  arrives  in  the  C- 
D  time  interval.  However,  it  is  actually  the  case  that  the  initial  arrival  is  in  the  A-B  window.  The 
low  noise  level  to  the  right  is  a  coincidence;  and  this  is  seen  in  a  later  section  by  examination  of 
the  similar  signal  from  the  event  Petit,  which  took  place  near  Tanana  and  was  also  recorded  at 
Palmyra  on  Small  Fry.  We  will  see  that,  although  the  signal  is  very  similar,  as  would  be  expected, 
the  amplitudes  after  the  window  corresponding  to  C-D  are  not  low. 
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DATA  ANALYSIS  TECHNIQUES 


To  compute  the  cross  correlations  of  the  signals  in  Figure  3,  we  first  align  them  in  time,  using  the 
delays  which  give  a  maximum  in  the  unfiltered  trace  cross-correlation  for  the  C-D  window  shown 
in  Figure  3.  This  window  is  chosen  because  it  appears  to  have  the  highest  S/N.  To  determine  the 
alignment  delays  we  note  that  the  signal  arrives  first  on  channel  2.  For  maximum  correlation, 
channel  1  must  be  moved  3.6  seconds  to  the  left,  and  channel  3, 3.3  seconds  to  the  left.  In  subse¬ 
quent  calculations  of  correlation,  we  calculate  the  correlation  at  these  fixed  offsets  for  both  the  A- 
B  and  C-D  windows,  and  for  all  periods. 

The  noise  mean-square  covariances,  needed  to  correct  the  correlation  estimates  as  discussed 
below,  were  calculated  for  each  channel  over  the  window  comprising  the  first  3  minutes  of  data  in 
Figure  3.  Cross-correlations  were  also  calculated  over  these  windows  and  were  found  to  be  ran¬ 
domly  positive  and  negative,  and  not  significantly  different  from  zero,  showing  that  there  was  no 
significant  signal  in  the  first  3  minutes  so  that  the  mean-square  covariance  estimates  were 
unbiased. 

The  correlations  are  calculated  using  the  above  fixed  delays  for  all  periods,  rather  than  adjusting 
the  delay  for  maximum  correlation  at  each  period,  because  analysis  suggests  that,  on  the  one 
hand,  the  hypothesis  of  a  fixed  delay  cannot  be  rejected  and  that,  on  the  other  hand,  at  some  peri¬ 
ods  the  S/N  is  poor  enough  that  using  the  maximum  correlation  as  the  signal  correlation  estimate 
would  result  in  an  positively  biased  estimate  for  correlation. 

On  the  other  hand,  in  the  presence  of  noise,  the  normalized  signal  cross-correlation  at  a  fixed 
delay  will  be  biased  low.  This  is  because  the  single-channel  covariances  which  are  multiplied 
together  (and  the  square  root  taken)  in  the  denominator  of  the  conventional  cross-correlation  esti¬ 
mator  are  biased  above  the  value  which  they  would  have  for  signal  alone  in  the  absence  of  noise. 
On  the  other  hand,  the  cross-variance  term  in  the  numerator  is  unbiased  by  the  noise.  This  is 
because  infrasound  noise  is  uncorrelated  with  both  the  signal  and  with  noise  on  other  channels. 

Assuming  for  convenience,  without  loss  of  generality  (for  unequal  noise  values  we  take  As  to  be 
the  square  root  of  the  product  of  the  power  of  the  two  channels),  that  the  noise  has  the  same  spec¬ 
trum  on  all  channels,  then  the  correlation,  p,  is: 


P  = 


where  Cy  is  the  cross  correlation  between  two  channels  and  where  As  is  the  autocorrelation.  The 
signal  cross-correlation  corrected  for  noise,  under  the  assumption  that  the  noise  is  stationary  and 
uncorrelated  with  the  signal,  is  given  by: 


Pc 


ciJ 


As  ~  An 


(3) 
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where  An  is  the  autocorrelation  before  the  signal.  Note  that  since  statistical  fluctuations  are 
expected  between  the  noise  before  the  signal  and  the  noise  during  the  signal,  it  is  possible  for  the 
noise-corrected  signal  correlation  estimate  to  be  greater  than  1 .0.  Also,  it  is  possible  that  the  cor¬ 
relation  estimate  may  be  undefined.  For  low  S/N,  it  may  sometimes  happen  that  the  pre-signal 
mean-square  noise  estimate  is  greater  than  the  estimate  for  mean-square  signal  plus  noise.  Thus, 
the  estimate  for  the  mean-square  noise  in  the  signal  window  may  be  negative. 

To  estimate  detection  thresholds,  and  standard  errors  of  azimuth  estimates  at  those  thresholds,  it  is 
necessary  to  be  able  to  predict  the  S/N  and  correlation  for  yields  lower  by  a  factor,  Ay  than  those 
events  that  provided  the  basic  data. 

We  need  estimates  of  the  total  signal  power,  the  correlated  signal  power,  and  the  uncorrelated  sig¬ 
nal  power. 

It  is  interesting  to  note  that  the  uncorrelated  signal  power  does  increase  the  S/N  of  the  beam  in  the 
time  domain.  If  P  is  the  number  of  channels  in  the  beam,  the  uncorrelated  power  is  reduced  in  the 
beam  by  a  factor  of  1/P  in  comparison  to  the  correlated  power.  So  this  additional  uncorrelated 
power  is  added  to  the  correlated  power  after  the  arrival  of  the  signal,  enhancing  the  S/N  over  that 
due  to  the  correlated  power  alone. 

We  see  that  the  correlated  signal  power  on  a  single  channel  of  a  channel  pair  is: 

As,  c  ~  Pc'(As~An)  (4) 

and  the  uncorrelated  signal  power  on  a  single  channel  of  a  channel  pair  is: 

As,uc  =  <5> 

To  examine  a  case  of  a  reduced  amplitude  source,  we  may  multiply  Asc  and  As  uc  by  the  same  suit¬ 
able  source  strength  factor,  Ay,  to  obtain  A's  c  and  A's  uc‘,  here  uc  and  c  refer  to  uncorrelated  and 
correlated,  respectively.  Then,  the  correlation  would  be  given  by  the  cross-spectral  power  divided 
by  the  total  power  on  either  (equal  power)  single  channel,  comprising  the  sum  of  the  correlated 
and  uncorrelated  signal  power  plus  the  noise  power.  Thus, 


P' 


A’ 


s,  c 


A’s,  c  +  A's,  uc  +  An 


(6) 


would  be  the  correlation  between  the  same  two  channels  for  the  same  two  signal  waveform 
shapes,  but  with  each  amplitude  changed  by  the  factor,  Ay  and  embedded  in  the  original  noise,  An. 
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For  application  to  determine  probability  of  detection  and  errors  of  azimuth  estimates,  using  theo¬ 
retical  correlation  estimates  for  reduced  amplitude  signals  in  the  original  noise,  we  also  need  the 
signal-to-noise  on  an  individual  channel,  (S/N)-P  Thus, 


<.S/N)'p  =  (A’S'C  +  A’SiUC)/A 


(7) 


COMPARISON  OF  TANANA-PALMYRA  DATA 
TO  DATA  FROM  MACK  AND  FLINN  (1971) 


Table  1  gives  the  raw  signal  correlation  and  the  correlation  corrected  for  noise  for  the  two  win¬ 
dows  seen  in  Figure  3  and  for  the  bandpass  limits  illustrated  in  Figure  4.  The  actual  calculations 
were  performed  using  SAC2000  and  following  the  formulas  given  in  the  previous  section. 


Cross-Correlation  of  Aligned  Traces 

1:2, 1:3, 2:3 

Tanana  at  Palmyra  (Small  Fry) 

B  (Hz) 

Observed  Signal  Correlation 

Corrected  for  Background  Noise 

T  (sec) 

Window  A-B 

Window  C-D 

Window  A-B 

Window  C-D 

.025-0.1 

-.02 

.85 

U 

1.18 

20 

.35  .29 

.94  .85 

4.3  U 

1.36  1.30 

0.05-0.2 

.72 

1.3 

1.2 

10 

.34  .46 

.88  .87 

.54  .96 

1.3  1.3 

0.1 -0.5 

.72 

.69 

.90 

1.04 

4 

.62  .84 

.59  .68 

.80  .85 

.81  .86 

0.2-0.5 

.64 

.96 

.81 

1.4 

3 

.70  .82 

.63  .68 

.83  .90 

.98  .82 

0.3-1. 0 

.58 

.84 

.67 

.96 

2 

.48  .69 

.22  .60 

.55  .72 

.26  .64 

0.5-2.0 

.38 

.71 

.44 

.77 

1 

.28  .40 

.02  .21 

.35  .44 

.03  .22 

Table  1.  Correlation  Analysis  of  Signals  in  Figure  3 


Since  the  signals  arrive  at  channels  1  and  3  within  0.3  seconds  of  each  other,  we  may  assume  that 
the  vector  between  them  was  aligned  nearly  parallel  to  the  wavefront  and  that  they  are,  therefore, 
expected  to  have  the  lowest  correlation,  as  shown  by  MF.  The  greatest  delay,  of  3.6  seconds,  is 
between  channels  1  and  2  and  so  we  expect  that  these  are  oriented  most  perpendicular  to  the 
wavefront  and  would,  therefore,  have  the  highest  correlation. 

Examination  of  Table  1  shows  that  these  expectations  are  bom  out.  Assuming  that  all  undefined 
correlation  values,  and  corrected  values  greater  than  1 .0  (which  are  long-period)  may  be  approxi- 
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mated  by  1.0,  and  combining  values  for  both  windows,  we  have,  for  interpolation  purposes,  fit  the 
correlation  data  with  the  expression  1  -exp(-aT)  where  T  is  the  period  in  seconds.  We  find,  parallel 

to  the  wavefront  a= 0.39  sec*1,  and  perpendicular  to  the  wavefront,  a=0.87  sec*1,  in  accordance 
with  the  observation  that  for  a  fixed  period  and  distance  between  elements  the  correlation  is 
higher  perpendicular  to  the  wavefront  than  parallel  to  the  wavefront.  (We  assume  that  even 
though  the  signals  arrive  at  almost  the  same  time  at  elements  1  and  3,  the  array  is  roughly  in  the 
shape  of  an  equilateral  triangle  and  that  the  spacing  is  approximately  1  km.) 

Using  these  formulas,  we  may  compare  in  Table  2  the  interpolated  results  with  the  results  pre¬ 
dicted  by  MF.  The  predicted  and  observed  values  seem  to  be  in  reasonable  agreement  so  we  are 
encouraged  to  continue  to  use  the  formulas  of  MF  in  order  to  estimate  signal  correlation  as  a  func¬ 
tion  of  array  element  separation  azimuth  and  distance,  and  as  a  function  of  signal  period. 


Comparison  of  Predicted  and  Observed  Infrasound  Correlation 
for  Sensor  Spacing  of  1  km 

Coherence  Prediction, 
Mack  and  Flinn  (1971) 
Ac=0.015  km/sec 

A6=5° 

Fitted  Correlation 
Observations 

Tanana  at  Palmyra 

Signal  Period  (sec) 

2 

1 

2 

1 

p  parallel  to  wavefront 

.87 

.53 

.54 

.32 

p  perpendicular  to  wavefront 

.96 

.84 

.82 

.58 

Table  2.  Twelve  (12)  observations  from  T=  1  to  20  seconds,  from  Table  1,  constrained  to  be  <=1, 
were  fitted  to  the  expression  l-exp(-aT).  a  was  found  to  be  0.39  parallel  to  the  wavefront,  and 
0.87  perpendicular.  Considering  the  scatter  in  the  small  amount  of  data  which  are  available  for 
analysis,  the  observations  seem  not  inconsistent  with  the  coherence  predictions  of  MF  although 
there  is  an  indication  that  the  observed  correlation  is  less  than  predicted,  which  would  suggest  the 
use  of  smaller  arrays  than  those  suggested  using  the  data  of  MF. 
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DETECTION  THEORY 


The  theory  required  to  estimate  the  errors  in  azimuth  as  a  function  of  array  aperture  and  S/N  was 
given  in  Blandford  (1997),  Shumway  et  al.  (1999)  and  is  reproduced  in  Appendix  I.  However, 
evaluation  of  an  array  design  or  aperture  should  be  evaluated  in  terms  of  detection  as  well  as  azi¬ 
muth  estimation.  To  evaluate  detection  from  a  quantitative  point  of  view,  we  need  a  detection 
theory. 

Because  of  the  frequent,  large,  uncorrelated  bursts  of  wind  noise  which  appear  on  elements  of 
infrasound  arrays,  the  standard  STA/LTA  detectors  conventionally  used  for  seismic  detection, 
cannot  be  used  for  infrasound  detection;  their  use  would  result  in  a  very  high  false  alarm  rate. 
Instead,  detectors  have  been  used  for  decades  which  rely  on  detecting  the  presence  of  the  same 
signal  on  different  channels  of  the  array;  that  is,  detectors  which  rely  on  the  detection  of  correla¬ 
tion  between  the  channels.  See,  for  example,  Jacobson  (1957),  Smart  and  Flinn  (1971),  McKissic 
(1996). 

The  “F  detector”  discussed  by  Smart  and  Flinn  (1971)  and  by  Blandford  (1974)  does,  in  fact,  rely 
on  correlation  between  sensors,  and  also  has  a  detection  theory  which  is  firmly  based  in  classical 
statistics. 

The  F  detector  detects  on  the  ratio  of  the  array  beam  power  to  the  mean-over-channels  residual 
power.  The  residual  power  is  computed  by  squaring  the  result  obtained  by  subtracting  the  beam 
from  the  aligned  signal  channel  and  then  averaging  over  channels.  The  ratio,  a  ratio  of  chi-square 
statistics,  is  an  F  statistic. 

When  a  signal  arrives,  the  ratio  of  chi-squares  increases  for  two  reasons.  First,  the  numerator 
increases  because  the  beam  signal  power  is  added  to  the  ongoing  beam  noise  power.  Second,  if 
the  signal  is  correlated,  then  when  the  beam  is  subtracted  from  each  individual  channel,  the  resid¬ 
ual  noise  is  the  same  as  before  the  signal  arrived,  and  the  numerator  remains  the  same  so  that  the 
ratio  increases. 

However,  if  an  uncorrelated  signal  or  burst  of  noise  arrives,  then  the  beam  will  increase  in  ampli¬ 
tude,  but  by  less  than  if  the  signal  were  coherent,  and,  in  the  numerator,  subtracting  the  beam 
increases  the  residual.  As  a  result,  the  numerator  and  denominator  increase  in  parallel  and  the 
ratio  remains  much  the  same. 

In  the  presence  of  signal,  the  ratio  is  distributed  as  non-central  F.  The  cumulative  integral  of  F 
from  0  to  q  with  Nj  and  N2  degrees  of  freedom  and  non-centrality  parameter  X  is  expressed  as: 

pF(q,N1,N2,X)  (8) 

where  N^BT  is  the  degrees  of  freedom  for  the  numerator,  where  B  is  the  bandwidth  (for  exam¬ 
ple  one  of  those  indicated  in  Table  1),  and  T  is  the  length  of  the  detection  window  over  which  the 
power  is  averaged.  In  this  report,  we  take  T=100  seconds,  which  Figure  3  shows  to  be  suitable. 
For  the  denominator,  N2=(P-1)2BT  where  P  is  the  number  of  elements  in  the  array. 
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The  non-centrality  parameter,  X  is  given  by  X=2BT(S/N)beam,  where  (S/N)beam  contains  the 
correlated  power  (the  true  signal,  S)  and  is  the  S/N  on  the  beam. 

For  detection  analysis,  it  is  necessary  to  fix  the  false  alarm  rate.  This  corresponds  to  the  probabil¬ 
ity  of  the  F  statistic  fluctuating,  in  the  absence  of  signal  (X=0),  and  exceeding  a  fixed  threshold, 
qT,  in  a  period  of  time.  We  choose  the  rate  to  be  one  false  alarm  per  day.  This  false  alarm  rate  per 
beam  would  result  in  approximately  10  false  alarms  per  day  for  a  200-meter  aperture  infrasound 
array  which  typically  has  10  noise-independent  beams.  This  number  of  false  alarms  is  expected 
to  result  in  very  few  false  two-station  network  events  per  day.  For  the  selected  detection  window 
of  T=100  seconds  (again,  seen  to  be  suitable  in  Figure  3),  this  corresponds  to  l-pF(qT,  Nl5  N2, 
0)=0.001 16.  Solving  for  qT  for  the  varying  values  of  B  and  for  P=4  gives  the  values  for  qj  found 
in  Table  3. 


Detection  and  Location  Parameters  for  CD  Infrasound  (1,2,3)  km  Aperture  Arrays  Derived  from 
Mack  and  Flinn  (1971)  Correlations,  and  Power  (S/N)p  from  Tanana  as  Observed  at  Palmyra, 

Reduced  by  a  Factor  of  Aj=l/25. 

Tanana  at  Palmyra  Signal  Parameters 

Theoretical  Mean 
Signal  Correlation,  p , 
for  Int’l  Monitoring 
Community  Arrays 
(all  numbers) 

Non-centrality  Parameter 
X=2BT*(S/N)beam 
for  Int’l  Monitoring 
Community  Arrays, 

Aj= 1/25 
(all  numbers) 

H 

B 

(Hz) 

qT 

P 

D 

S/Np 

1km 

2km 

3km 

1km 

2km 

3km 

20 

.075 

3.12 

0.54 

1.0 

1.2 

.998 

.995 

.99 

1.62 

1.62 

1.62 

10 

.15 

2.33 

0.62 

1.0 

1.6 

.995 

.98 

.96 

14.9 

14.6 

14.3 

D 

1.70 

0.69 

.87 

3.8 

.97 

.89 

.77 

55.2 

48.8 

42.0 

3 

.3 

1.84 

0.74 

.89 

4.9 

.95 

.81 

.63 

41.6 

34.8 

26.5 

2 

MM 

1.50 

0.57 

.63 

8.1 

.89 

.63 

.39 

130 

87.2 

50.4 

1 

1.5 

1.32 

0.33 

.38 

13.2 

.63 

.22 

.065 

262 

77.5 

21.6 

Table  3.  Detection  and  Location  Parameters  for  CD  Infrasound 


Tc  is  the  center  period  for  the  detection  filter  passband  with  bandwidth  B  and  time  duration  T=  1 00 
sec.  qT  is  the  F  detection  threshold  for  one  false  alarm  per  day.  p  is  the  average  observed  corre¬ 
lation,  which  is  estimated  as  the  average  over  the  data  in  columns  2-3  of  Table  1  for  the  observed 
correlation  between  two  traces. 
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(S/N)p  is  the  power  signal-to-noise  ratio  on  a  single  channel;  it  may  be  derived  starting  with  the 
relation,  for  the  observed  correlation  between  two  traces 

_  _  _ con^ejatedsigncd _ 

uncorrelatedsignal  +  noise  +  correlatedsignal 


where  we  may  take  the  correlated  signal  as  equal  to  the  signal  power  S  times  the  correlation  cor¬ 
rected  for  background  noise,  pc .  (We  may  estimate  pc  as  the  average  over  columns  4-5  in 
Table  1,  with  1.0  replacing  all  values  greater  than  1.)  Then,  the  uncorrelated  signal  may  be  deter¬ 
mined  as  (l-pc)  times  S. 

Using  these  relations,  one  may  obtain  the  formula, 


(10) 


This  formula  is  used  in  Table  3  to  produce  column  6  from  columns  4  and  5.  The  approximation 
here  is  that  we  average  over  relative  propagation  azimuths  and  assume  that  all  three  sensor  spac- 
ings  are  the  same  for  the  observing  array  at  Palmyra. 

( S/N)P  in  Table  3  is  one  of  the  critical  pieces  of  information  determined  by  this  analysis  of 

Tanana  as  recorded  at  Palmyra.  It  gives  an  estimate  of  the  single-sensor  power  S/N  as  a  function 
of  frequency  for  a  source,  distance,  and  system  of  interest. 

It  should  be  noted  that,  as  mentioned  previously,  due  to  the  loss  of  signal  amplitude  created  by 
300-meter  pipe  arrays,  it  could  be  that  this  value  of  S/N  at  1  Hz  is  underestimated  compared  to 
that  which  would  be  possible  for  a  more  ideal  system.  In  addition,  it  is  possible  that  the  space-fill¬ 
ing,  100m  aperture  pipe  arrays  currently  being  considered  for  the  International  Monitoring  Com¬ 
munity  arrays  may  improve  S/N  over  that  found  here. 

However,  going  ahead  and  using  these  results,  together  with  the  theoretical  relations  for  correla¬ 
tion  as  a  function  of  period  and  element  spacing  from  MF  (also  verified  by  the  Tanana-Palmyra 
data),  we  will  calculate  S/N  on  the  beam,  detection  probability,  and  azimuth  estimation  precision 
for  arrays  of  different  designs  looking  at  the  Tanana  signal  scaled  to  simulate  a  smaller  yield. 

The  columns  headed  “Theoretical  Mean  Signal  Correlation”  in  Table  3  are  determined  using  the 
parameters  Ac  =0.0 1 5  km/sec  and  A0  =5°  and  the  relation  of  MF  to  determine  the  mean  correla¬ 
tion  for  signals  at  infinite  S/N  at  International  Monitoring  Community  arrays  with  apertures  of  1, 
2,  and  3  km. 

Note  that  the  theoretical  values  for  a  1  km  array  in  Table  3,  column  7,  are  larger  than  are  the  cor¬ 
rected  observed  values  in  column  5,  pc .  This  is  appropriate  because  three  of  the  six  intersensor 
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distances  in  a  4-element  International  Monitoring  Community  array  (an  equilateral  triangle,  1  km 
on  a  side,  plus  a  central  element)  of  this  size  are  0.43  km;  whereas  the  3-element  Small  Fry  array 
had  only  three  spacings,  each  likely  was  approximately  1  km.  Also,  we  have  noted  above,  in  dis¬ 
cussion  of  Table  2,  that  the  observed  correlations  appear  to  be  somewhat  lower  than  predicted. 

For  detection  analysis  purposes,  it  is  necessary  to  calculate  the  non-centrality  parameter,  X .  The 
new  element  in  this  calculation  is  given  by  the  (S/N)  on  the  beam.  To  be  consistent  with  the 
classical  development  of  the  F  detector,  Shumway  (1983),  which  assumes  that  the  signal  is  per¬ 
fectly  correlated,  we  assume  that  only  the  correlated  signal  power  will  count  as  S.  Thus,  counting 
the  uncorrelated  S  as  noise,  (S/N)t,eam  is  given  for  a  4-element  array  by 


:s/N)beam  = 


p  -AjJJVN)p 
\-[{\-p)-Af'{S/N)p+\] 


(11) 


where  Aj-is  the  factor  giving  the  reduction  of  signal  power;  for  Table  3,  Af  -  1/25,  p  is  the  average 
corrected  signal  correlation  for  the  array  (columns  7-9  in  Table  3).  This  is  equivalent  to  (l/25)th 
of  the  yield  or  0.104  kt.  We  test  capability  for  this  low  yield  in  order  to  evaluate  in  a  yield  region 
where  there  are  significant  operational  differences  between  different  arrays  sizes.  For  substan¬ 
tially  larger  yields,  all  arrays  would  perform  adequately  at  this  distance,  and  for  much  smaller 
yields,  no  arrays  would  perform  adequately. 
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DETECTION  AND  AZIMUTH  ESTIMATION  RESULTS 


Using  the  parameters  determined  in  Table  3,  together  with  the  detection  theory  outlined  in  the  pre¬ 
vious  section,  and  with  the  azimuth  estimation  theory  in  Blandford  (1997)  (see  also  Appendix  I), 
we  obtain  the  results  in  Table  4. 


Detection  and  Location  Estimates  for  CD  Infrasound  (1,2,3)  km  Aperture  Arrays 

Using  Parameters  Derived  from  Mack  and  Flinn  (1971),  and  the  Event  Tanana  as 
Observed  at  Palmyra.  Estimates  for  the  Event  Tanana  as  Observed  at  Palmyra  with 

Amplitude  Reduced  by  a  Factor  of  5. 

Tc(sec) 

Beam  Visual 
Amplitude  (S/N) 
(rms/rms) 

Probability  of  Detection 
With  F  using  Simple 
Beam 

PD 

Azimuth  Estimate  with 
Optimum  Weighted 
Beam 

Standard  Error,  Og  0 

1  km 

2km 

3  km 

1  km 

2  km 

3  km 

1km 

2km 

3  km 

20 

1.09 

1.09 

1.09 

.006 

.006 

.006 

20.1 

10.1 

6.8 

10 

1.12 

1.12 

1.12 

.012 

.012 

.012 

6.10 

3.16 

2.23 

4 

1.26 

1.26 

1.24 

.34 

.31 

.25 

1.04 

.69 

.62 

3 

1.33 

1.31 

1.29 

.39 

.32 

.22 

.89 

.68 

.63 

2 

1.49 

1.44 

1.37 

.998 

.977 

.78 

.42 

.36 

.36 

1 

1.66 

1.46 

1.36 

1.0 

.98 

.41 

.21 

.31 

.66 

Table  4.  Detection  and  Location  Estimates  for  CD  Infrasound 


In  Table  4,  Tc  is  the  center  period  for  the  detection  filter  passband.  Detection  and  azimuth  esti¬ 
mates  are  for  a  time  window  of  100  seconds.  The  “Visual”  amplitude  in  columns  2-4  is  computed 
by  dividing  the  estimate  for  the  beam  noise  power  in  front  of  the  signal  into  the  sum  of  the  beam 
noise  plus  correlated  signal  power  plus  1/4  the  uncorrelated  signal  power  (for  a  4-element  array), 
and  then  taking  the  square  root. 

The  probability  of  detection  is  computed  using  the  Splus  statistical  package  pF  subroutine  which 
is  suitable  for  a  correlation  F  detector.  A  false  alarm  rate  of  one  per  day  is  specified  along  with  a 
signal  with  l/25th  the  power  (yield~0.104kt)  of  Tanana  at  Palmyra.  The  azimuth  estimate  was 
calculated  from  the  formulas  in  Appendix  I  using  a  Matlab  analysis  package  routine,  fsnewth.m, 
written  by  the  author. 

Calculation  of  azimuth  error,  of  course,  presupposes  detection  and  also  presupposes  that  subse¬ 
quent  analysis  has  selected  the  proper  time  window  over  which  to  estimate  the  azimuth.  Thus,  if 
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detection  probability  is  low,  the  significance  of  the  azimuth  estimate  is  questionable.  Examination 
of  Table  3  appears  to  suggest  that  one  should  construct  1  km  aperture  arrays  and  detect  at  1  Hz. 
Note  that  the  best  detection  and  azimuth  estimate  for  this  signal  is  found  at  1 -second  period  for  a 
1  km  aperture  array. 

For  1  and  2  km  aperture  arrays,  the  best  detection  is  found  at  1 -second  period.  For  3  km  aperture 
arrays,  the  best  detection  is  found  at  2-second  period. 

Using  the  best  period  for  detection  for  each  aperture  array,  we  may  find  the  yield  which  gives  90% 
probability  of  detection  for  a  fixed  false  alarm  rate  of  one  per  day.  The  result  is  0.032, 0.078,  and 
0.124  kt  for  1, 2,  and  3  km  aperture  arrays,  respectively.  Thus,  a  1  km  aperture  array  will  detect 
1/2.4  the  yield  as  a  2  km  aperture  array,  and  1/3.9  the  yield  as  a  3  km  aperture  array. 

The  theory  and  expressions  for  the  azimuthal  error  when  simple  beamforming  is  used  instead  of 
optimal  beamforming  are  also  given  on  the  last  page  of  the  Appendix.  Perhaps  surprisingly,  both 
the  theory  and  the  expressions  are  somewhat  more  complex  for  simple  beamforming.  However, 
again  the  expression  for  the  azimuthal  error  has  been  computed  for  simple  beamforming  using  the 
Matlab  subroutine,  bfsnewth.m  (see  Appendix),  also  written  by  the  author. 

In  the  last  row  of  Table  4,  for  1 -second  period,  the  azimuthal  error  in  degrees  for  (1, 2, 3)  km 
aperture  for  simple  beamforming  would  be  (0.33, 1 .00, 1 .32)  instead  of  the  (0.21 5, 0.3 1 , 0.66) 
seen  in  the  table,  which  is  for  optimum  processing.  For  2-second  period,  the  values  are  (0.44, 
0.51, 0.88)  as  compared  to  (0.42, 0.36, 0.36),  and  for  longer  periods,  the  differences  are  30%  or 
less. 

Thus,  for  the  optimum  1  km  aperture  array,  there  is  only  a  1.5  factor  (0.33/0.215)  advantage  in 
azimuth  precision  to  using  the  optimum  azimuth  estimation  process  in  preference  to  fk  beam¬ 
forming  at  the  optimum  detection  period  (1  sec),  while,  for  2  km  and  3  km  aperture  arrays,  the 
advantages  are  factors  of  3.2  (at  1  sec)  and  2.4  (at  2  sec),  respectively. 

Thus,  the  1  km  aperture  array  offers  an  advantage  in  processing  in  that  there  is  a  lesser  require¬ 
ment  for  processing  that  is  more  sophisticated  than  simple  beamforming  to  obtain  the  full  benefit 
of  the  data  than  there  is  for  the  2  km  and  3  km  arrays. 

Note  also  that  we  have  used  as  a  detector  the  F  statistic  on  a  simple  beam.  It  may  be  that  the 
detection  performance  of  the  larger  arrays  could  be  improved  by  an  optimum  detector  which  took 
account  of  the  loss  of  signal  coherence.  In  a  qualitative  sense,  such  a  detector  would,  presumably, 
detect  to  some  degree  on  the  sum  of  the  power  at  each  sensor.  Such  a  detector  would,  however, 
likely  be  very  vulnerable  to  the  bursts  of  incoherent  noise  which  are  characteristic  of  inffasound 
data.  Thus,  again,  we  are  brought  back  to  using  the  smaller  1  km  aperture  array. 

Of  course,  we  must  emphasize  that  all  these  conclusions  rest  on  the  correlation  predicted  by  Mack 
and  Flinn  (1971),  which  was  confirmed  by  Tanana,  and  on  the  observed  (S/N)  at  a  single  sensor  as 
a  function  of  frequency  for  the  event  Tanana  as  seen  at  Palmyra.  Further  data  is  sorely  needed. 
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ADDITIONAL  CORRELATION  DATA 


The  data  for  the  event  Tanana  as  recorded  at  Small  Fry  is  the  best  data  which  I  have  been  able  to 
acquire  for  the  purpose  of  testing  whether  or  not  the  MF  data  could  be  extended  to  shorter  periods 
and  shorter  distances  between  sensors. 

However,  there  are  some  further  data  which  are  consistent  with  the  Tanana-Small  Fry  data.  In 
Figure  5,  we  see  the  event  Petit,  also  recorded  at  Tanana.  Analysis  of  the  same  signal  windows  as 
seen  in  Figure  3  gave  similar  results  but  with  larger  variation. 

Figure  6  shows  Petit  as  recorded  at  a  10-km  aperture  array  at  Oahu,  a  distance  of  approximately 
2200  km.  Unfortunately,  the  archived  signal  did  not  include  a  portion  of  uncorrelated  noise 
before  the  signal  arrival  so  that  it  was  not  possible  to  determine  the  true  correlation  values  for 
these  signals,  only  a  lower  value.  These  values,  together  with  the  theoretical  MF  values,  are  seen 
in  Table  5  below. 


Comparison  of  Observed  and  Theoretical  Infrasound  Signal  Correlations 

Period  (sec) 

20 

10 

n 

3 

2 

1 

10  km  spacing,  Petit,  2.2  kt,  at  Oahu  (-2200  km) 

Perpendicular 
to  wavefront, 
elements  1-2 

Observed  (raw) 
(:13:50-:15:30  window) 

>.91 

>.56 

>.29 

>.35 

>.13 

>-.09 

Theoretical 

.96 

.84 

.24 

.05 

-0 

~o 

Parallel 
to  wavefront, 
elements  1-3 

Observed  (raw) 
(:13:50-:15:30  window) 

>.68 

>.41 

>-.15 

>-.28 

>-.25 

>-.09 

Theoretical 

.87 

.53 

.22 

.03 

-0 

-0 

Table  5.  Signal  Correlation  Comparison 


If  the  negative  values  of  correlation  parallel  to  the  wavefront  for  4, 3, 2,  and  1 -second  periods  in 
Table  5  can  be  taken  to  indicate  that  the  standard  errors  of  the  Oahu  correlation  estimates  are 
-0.25,  then  the  theoretical  (Mack  and  Flinn,  1971),  and  observed  data  seem  not  inconsistent, 
although  there  is  a  suggestion  that,  for  estimates  perpendicular  to  the  wavefront,  the  observed  cor¬ 
relations  are  greater  than  the  theoretical. 

Figure  7  shows  the  White  Sands  19  November  1997  10-ton  conventional  explosion  as  seen  at  the 
TXAR  infrasound  array.  In  general,  the  S/N  is  fair.  Figure  8  shows  the  01  channel  filtered 
through  several  passbands,  and  we  see  that  a  signal  appears  to  be  present  between  20  and  1/2 
second  period.  ! 
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The  White  Sands  event  is  approximately  450  km  to  the  northwest,  such  that  only  the  1.1  km  inter¬ 
val  between  elements  01  and  02  is  roughly  perpendicular  to  the  wavefront.  The  intervals  (01  - 
03),  (01  -  09),  and  (03  -  09)  are  roughly  parallel  to  the  wavefront  and  are  0.8, 1 .5,  and  2.0  kilome¬ 
ters  in  length,  respectively. 

Table  6  gives  correlations  for  these  distances,  calculated  in  the  same  way  as  for  those  in  Table  1 
and  for  a  200-second  window  surrounding  18:27  as  seen  in  Figures  7  and  8. 


Cross-Correlation,  pc , corrected  for  noise, 

of  Aligned  Traces 

White  Sands,  10  tons,  at  TXAR 

01 :02  perpendicular  to  wavefront,  others,  parallel 

B  (Hz) 

perpendicular 
intervals.  1  km 
01:02 

parallel  ! 

interval=0.8  km 
01:03 

T  (sec) 

parallel 

interval=1.5  km 
01:09 

parallel 

interval=2.0  km 
03:09 

2.0-4.0 

.75 

54 

.5 

.79 

.46 

1. 0-4.0 

.97 

.83 

0.7 

.89 

.88 

0.5-2.0 

.97 

.97 

1.5 

.97 

.98 

Table  6.  Correlation  Analysis  of  Signals  in  Figures  7  and  8  for  Four  Different  Sensor  Pairs. 

Periods  given  in  Table  6  are  the  observed  dominant  period  of  the  signal  as  measured  after  filtering. 
It  is  apparent  that  only  correlations  greater  than  0.97  would  be  observed  for  periods  greater  than 
1 .5  seconds. 

At  periods  shorter  than  0.5  seconds,  the  S/N  was  poor,  as  can  be  seen  in  Figure  8.  This  is  for  a 
small  (10  tons,  0.005  kt)  shot,  short  distances,  and  with  a  hose  array  with  a  diameter  of  only  20 
meters.  Thus,  it  seems  unlikely  that  higher  frequencies  will  be  useful,  given  the  capability  at 
lower  frequencies,  for  detection  or  azimuth  estimation  of  larger  explosions  at  larger  distances.  It 
may  possibly  be  that  the  higher  frequencies  would  be  useful  for  discrimination  if  they  are  gener¬ 
ated  in  greater  proportion  by  nonexplosive  sources. 
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It  is  important  to  emphasize  that  these  high  correlation  values  would  be  lower  if  not  corrected  for 
noise.  For  example,  the  0.98  value  for  the  1 .5-second  period  for  the  03:09  element  pair  separated 
by  2  km  parallel  to  the  wavefront  would  be  0.89  without  correction;  still  quite  high,  of  course. 

Comparison  of  the  squared  values  of  the  1.5  second  correlations  with  Figure  1  shows  that  such 
high  correlations  at  distances  between  sensors  of  1-2  km  would  only  be  obtained,  with  MF  param¬ 
eters,  at  periods  of  1 0  seconds.  This  implies,  roughly,  that  the  azimuth  and  phase  velocity  vari¬ 
ance  determined  by  MF  would  have  to  be  reduced  by  about  a  factor  of  5  in  order  to  be  applicable 
to  these  regional  signals. 

One  interpretation  of  this  result  might  be  that  a  one-hop  signal,  present  at  closer  distances,  does 
not  have  the  azimuthal  and  velocity  dispersion  that  a  multi-hop  signal  has. 

If  true,  this  result  emphasizes  that  results  appropriate  for  1000-2000  km  detection  and  azimuth 
estimation  are  unlikely  to  be  obtained  by  analysis  of  records  from  500  km  distances.  Further 
work  is  required. 
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SUMMARY  AND  CONCLUSIONS 


It  seems  likely  that  the  optimum  design  for  infrasound  arrays  for  detection  and  location  of  1  kt 
nuclear  explosions  at  1 000-2000  km  should  be  1  km  aperture  arrays,  and  that  detection  processing 
should  be  concentrated  in  the  neighborhood  of  1  Hz.  This  conclusion  is,  if  anything,  likely  to  be 
strengthened  by  future  analyses  since  the  1  Hz  signal  recorded  at  Small  Fry  was  likely  degraded 
in  S/N  by  the  large  aperture  pipe  array. 

There  remains  a  paucity  of  data  for  a  few  kiloton  explosions  recorded  at  1000-2000  km.  There 
are  some  additional  White  Sands  4  kt  conventional  explosion  data  recorded  at  great  distances,  and 
substantial  efforts  to  recover  that  data  would  be  worthwhile.  At  AFTAC,  Dean  Clauter  (personal 
communication)  is  striving  to  recover  from  the  archives  further  nuclear  explosion  infrasound  data 
recorded  at  small  Pacific  Island  arrays.  These  data  should  be  analyzed  to  ensure  that  the  above 
conclusion  is  reliable. 

It  would  appear  that  smaller  explosions  at  closer  distances  will  also  have  optimum  detection  at 
1  Hz,  and  that  the  correlation  will  be  excellent  at  1  Hz  for  a  1  km  aperture  array,  so  that  location 
accuracy  should  be  excellent.  It  might  be  that  a  higher  aperture  array  would  give  better  location 
for  these  closer  distances,  but  this  is  not  likely  to  be  a  serious  International  Monitoring  Commu¬ 
nity  problem. 
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Figure  2.  Map  of  Palmyra  Island,  5.88°N,  162.08°W  which  tats  the  site  of  the  Small  Fry  infrasound  array  which  recorded  the  signals 
of  the  event  Tanana  on  25  May  1962.  Tanana  was  detonated  25  km  south  of  Kiritimati  (Christmas)  Island  at  1.65°N,  157.28°W,  a  dis¬ 
tance  of  708  km  to  the  SE  of  Palmyra.  The  array  indicated  on. the  map  is  not  the  1962  3-element  Small  Fry  array,  whose  coordinates  are 
at  present  unknown,  but  is  instead  a  larger  array  installed  in  1 95 1 . 


Tanana,  2.6  kt,1.65N,  157.28W,  16:08:52,  25  May  1962  at  Palmyra  (Small  Fry) 


00:07:64.2  :10:00 


Signal  Windows:  1.A-B;  2,  C-D 


Figure  3.  Unfiltered  waveforms  of  Tanana  as  recorded  at  Palmyra.  The  first  arrival  is  shortly 
after  the  vertical  line  “A”.  The  two  signal  windows,  of  100  seconds  each,  are  the  A-B  and  C-D 
intervals.  A  1 80-second  noise  window  before  “A”  was  used  to  estimate  noise. 


23 


Tanana,  2.6  kt,1.65N,  157.28W,  16:08:52,  25  May  1962  at  Palmyra  (Small  Fry) 


23:57:29  00:00:00  :02:30  :05:00  :07:30  :10:00 


Time  (hr  mi  msec)  -  ~1 6:43:00 


Figure  4.  Channel  3  of  Tanana  as  recorded  at  Palmyra,  filtered  through  passbands  with  center 
periods  as  indicated.  The  S/N  is  roughly  equal  over  the  passbands;  however,  as  discussed  in  the 
text,  the  shorter  periods  have  the  advantage  of  greater  bandwidth  for  detection  and  azimuth  esti¬ 
mation  so  long  as  signal  correlation  is  adequate  for  the  array  sensor  spacing. 
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Pool,  19  June  62/Tan  25  May,  62 


25 


Figure  5.  Petit,  2.2  kt,  (code-name  Pool)  and  Tatiana,  2.6  kt,  (Tan)  as; recorded  at  the  array  Small  Fry  on  Palmyra.  The  figure  shows 
that  Tanana  had  a  better  S/N  ratio  and  that  the  sharp  reduction  in  noise  amplitude  after  the  second  arrival  for  Tanana  was  a  coincidence 
which  did  not  occur  for  Petit.  The  poorer  S/N  for  Petit  made  it  possible  only  to  show  that  Petit  was  not  inconsistent  with  the  results  for 
Tanana;  Petit  could  not  substantially  improve  the  correlation,  and  comparative  S/N  as  a  function  of  period,  estimated  with  Tanana. 


Petit,  -15:00, 1 9  June,  1 962  at  Oahu  (Grape  Sugar),  -2200  km 


17:10:00  :12:30  :15:00  :17:30  :20:00 

Time  (hr:min:sec) 


Figure  6.  Petit,  2.2  kt,  as  recorded  at  the  Grape  Sugar  array  on  Oahu,  approximately  2200  km 
distant.  Analysis  shows  substantial  correlation  at  the  beginning  of  the  window,  which  was  digi¬ 
tized  from  an  archival  reproduction,  so  that  no  noise  window  is  available  with  which  to  calculate 
corrected  correlations.  The  typical  dimension  of  the  array  is  10  km  and  delays  are  correspond¬ 
ingly  as  great  as  30  seconds. 


White  Sands  Explosion,  19  Nov  1997,  atTXIAR,  ~450  km 
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1  B:22:00  :24:00  :26:00  :28:00  :30:00 

Time  (hr:min:sec) 


Figure  7.  The  White  Sands,  New  Mexico,  simultaneous  conventional  explosion  of  0.005  kt  on 
1 9  November  1 997  as  recorded  at  the  TXAR  infrasound  array  at  a  distance  of  approximately 
450  km.  Note  that  channel  1  has  lower  noise  than  the  other  channels,  while  channel  2  has  higher 
noise.  Unequal  conditions  of  this  sort  prevail  at  many  infrasound  arrays,  and  detection  and  azi¬ 
muth  estimation  techniques  should  be  developed  to  allow  for  these  differences. 
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Whit©  Sands  Explosion,  19  Nov  1 997.TXI01 , 
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Figure  8.  Channel  1  of  TXAR  from  the  White  Sands  event,  the  lowest  noise  channel  shown  in 
Figure  7,  filtered  through  passbands  with  center  periods  as  indicated.  The  highest  S/N  is  for  peri¬ 
ods  of  2  and  1  seconds  and,  as  discussed  in  the  text,  the  correlation  is  over  0.9  for  both  periods 
(much  higher  for  this  “regional”  signal  than  would  be  predicted  by  the  parameters  estimated  by 
Mack  and  Flinn  (1971)  for  longer  paths).  Thus,  we  expect  that  the  best  detection  and  azimuth 
estimation  would  occur  at  1  second  since,  at  that  period,  greater  bandwidth  is  available. 
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APPENDIX:  Detection  and  Estimation  for  Propagating  Stochastic  Signals 

Robert  H.  Shumway 
Division,  of  Statistics 
University  of  California,  Davis 
October  1,  1996 


1.  The  Stochastic  Signal  Model 

In  general  a  model  for  a  collection  of  signals,  observed  at  an  array  of  N  sensors  whose 
response  is  denoted  by  Vj(t),j  =  1,  •  •  • ,  N,  t  =  0, 1, . . . ,  T  —  1  is 

yj(t)  =  sj{t-Tj(e))+vj(t)  (1) 

where  s;(t)  are  random  signals,  assumed  to  be  different  on  each  channel,  and  Tj(9)  are 
time  delays  induced  by  a  propagation  pattern  indexed  by  the  wavenumber  coordinates 
0  _  (0lie2)>  which  are  nonlinearly  related  to  the  velocity  and  azimuth  of  a  propagating 
plane  wave.  We  assume  that  the  signals  are  stationary  processes  and  denote  the  N  x  N 
spectral  matrix  by  5  <  v  <  .5.  The  noise  processes  are  often  assumed  to  be 

independent  and  identically  distributed  across  the  array  but  have  a  stationary  correlation 
structure  over  time;  we  denote  the  common  spectral  density  by  P«( v).  with  v  denoting 
frequency  in  cycles  per  unit  time.  The  time  delays  are  expressed  m  terms  of  the  physical 
locations  rj  =  (rji,rj2)'  of  the  sensors  as 

Tj(0)  =  f  (2) 

over  a  set  of  frequencies  v  where  the  assumption  (2)  can  be  made.  We  may  also  consider 
a  version  of  (1)  which  assumes  a  common  stochastic  signal,  say  s(i),  on  all  channels,  i.e. 

yj(t)=s(t-Tj(e))+vj(t),  (3) 

and  we  refer  to  this  as  the  perfect  correlation  model. 

It  is  conventional  to  consider  the  above  model  in  the  frequency  at  a  collection  of  L 
frequencies,  say  i/i, . . . ,  vL  over  which  the  signal  and  noise  spectral  matrices  are  const  an  , 
say^ at  f3  and  P„Ijv  where  IN  denotes  the  N  x  N  identity  matrix.  Talcing  discrete  Fourier 

transforms  yields  the  approximation 

Yjt  =  Aj(9)Sjt  +  Vje  (4) 

over  a  set  of  frequencies  t  —  1 , .  -  • ,  L,  where 

Aj(9 )  =  exp{27nVTj} 

=  exp{2irirj-^}.  (3) 
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This  means  that  we  may  write  a  vector  form  of  (4)  in  the  frequency  domain  as 

Ye  =  G(9)St  +  Vt,  (6) 

where  Y*,  S*  and  V*  are  vector  transforms  of  the  observed  data,  signal  and  noise  respec¬ 
tively  and 

G(9)  =  diag{Ai(6),  A2(9), . . . ,  Aw(0)}  (7) 

is  an  N  x  N  matrix  with  the  Aj{9)'s  down  the  diagonal.  The  spectral  matrix  of  the 
observed  vector  Y  is  clearly 


/,(«)  =  G(9)f,G"(6)  +  /„/„  (8) 

where  C*  denotes  the  complex  conjugate  transpose  of  the  matrix  C.  The  mean  value  of 
the  vector  Y  t  is  zero  since  the  signal  and  noise  processes  axe  assumed  to  have  zero  means. 

A  special  case  of  interest  is  the  perfectly  correlated  signal  where  we  assume  that  or 
equivalently,  that  the  model  (3)  is 


Yj,  =  Aj(6)S,  +  Vjt, 

(9) 

which  is  stacked  in  the  form 

Y,  =  g(9)St  +  V, 

(10) 

where 

g  (t))  =  (A1(e),A2(e),...,AN(6))' 

(11) 

now  becomes  an  IV  x  1  vector  and  the  spectral  matrix  becomes 

/„(«)  =  P,  g(«)g*(«)  +  Pr.lt*- 

(12) 

It  should  be  noted  that  the  model  discussed  in  this  note  differs  from  the  usual  case 
where  we  regard  the  signal  as  being  fixed  and  unknown,  but  identical  between  sensors.  In 
that  case,  the  model  looks  exactly  like  Equation  (9),  with  the  signal  assumed  to  be  fixed 
and  unknown.  Hence  the  vector  Y  in  this  case  will  have  mean  g(9)S  and  spectral  matrix 
PnlN-  The  theory  for  various  proposed  estimators  in  this  case  has  been  covered  in  Hinich 
and  Shaman  (1972),  Hinich  (1981),  Wu  (1982),  Shumway  (1982)  or  Brillinger  (1985). 

For  the  stochastic  signal  case,  we  consider  estimation  of  the  signal  and  its  mean  square 
error  in  the  next  section  and  then  move  to  sections  on  maximum  likelihood  detection  and 
estimation  in  the  following  sections. 

2.  Signal  Estimation 

For  the  stochastic  signal  case,  it  is  clear  that, in  the  frequency  domain,  the  signal 
estimation  problem  can  be  solved  under  either  the  linearity  or  Gaussian  assumptions  by 
computing  the  conditional  expectation  of  the  signal  given  the  data  under  the  Gaussian  as¬ 
sumption.  We  consider  the  two  cases  corresponding  to  models  assuming  perfect  correlation 
and  more  general  correlation  structures. 
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2.1  Signal  Estimation:  The  Perfectly  Correlated  Signal 

The  perfectly  correlated  case  has  also  been  considered  by  Harris  (1990).  Suppose  we 
consider  the  model  given  by  (9)  under  the  assumptions  on  the  noise  and  spectral  matrices 
summarized  in  (11).  Using  the  Gaussian  assumptions  leads  to  estimating  St  by 


St(d)  =  £(S«|Y*) 

=  P,  g*(9)  (p, g*(0)g(0)  +  P„/jv)  Y, 

=  (g'WgW  +  £)  s'(«)y< 

_  rwY < 

N  +  r  ' 

where  g*(0)g(0)  =  N  from  (5)  and  (11)  and 

r  PM 
~  PM 


(13) 

(14) 


is  the  inverse  of  the  signal  to  noise  ratio.  The  result  exhibits  the  signal  estimator  as  the 
beam  g*(0)Y  adjusted  by  a  multiplier  that  depends  on  the  number  of  elements  N  and  the 
noise  to  signal  ratio  r.  The  mean  square  error  of  the  signal  estimator  reduces  to 


C T2{9 )  =  Pa~  P»g*(0)  ^Psg(^)  g*(0)  +  PnlN^j  gWPs 
=  P«(  g*(0)g(*)  +  7T 


Pn 

N  +  r 


(15) 


It  is  clear  that  the  estimated  signal  is  basically  the  beamforxned  estimator  g,*(0)Yt  and 
that  the  variance  goes  down  by  a  factor  that  is  weighted  by  the  number  of  sensors  plus 
the  inverse  of  the  signal  to  noise  ratio. 


2.2  Signal  Estimation:  The  General  Correlated  Signal 

For  the  general  stochastic  signal  case,  it  is  convenient  to  assume  that  the  transforms  are 
complex  Gaussian  and  again  use  the  fact  tnat  the  conditional  expectation  S(0)  £"(S^  j  ) 
has  the  smallest  mean  square  error.  For  the  general  model  (5),  we  obtain 


S*(0)  =  faG*(d)  (G(e)faG*(6)  +  PnlN^J  Y, 
=  (G*(6)G(9)  +  Pn/71)  G*(9)  Yt, 


(16) 
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using  a  standard  matrix  identity  (see,  for  example,  Jazwinski,  1970).  The  estimator  in¬ 
volves  computing  the  delayed  quantity  Aj(9)Yjt  on  each  sensor  and  then  adjusting  by 
multiplying  by  the  adjustment  matrix  involving  the  spectral  matrices  of  the  signal  and 
noise.  The  mean  square  covariance  matrix  of  the  estimator  is  given  by 

2(0)  =  /.-  f,G'($)(G(O)f,G'(0)  +  PnlA  G(B)f. 

=  pjG'(6)G(f>)  +  Pn/7')  ,  (17) 

using  another  identity  from  Jazswinski  (1970).  We  note  that  for  Aj(9)  as  defined  in  (3) 
and  (4),  we  have  the  simplification 


G*(9)G(9)  =  IN,  (18) 

so  that  the  multilplying  matrices  in  (16)  and  (17)  do  not  depend  on  9  and  we  may  write 


S  t{9)  =  CG*Ye 

(19) 

and 

S  =  PnC , 

(20) 

where 

C  =  (lN  +  PJ71  J  , 

(21) 

for  use  in  later  equations.  We  notice  that  the  optimal  estimator  is  no 
but  is  essentially  a  weighted  beam  of  the  form 

longer  the  beam, 

N 

Sji(9)  -  ^  cik  exp{2nir'k9}Ykl 

k= 1 

(22) 

with  weights  proportional  to  the  elements  of  C  defined  by  (21). 

3.  Maximum  Likelihood  Estimation 

In  the  stochastic  signal  case,  we  regard  the  wavenumber  vector  9  = 
nonlinear  functions  velocity  and  azimuth,  say 

=  (#i>  62)'  and  the 

V 

V 

=  m 

(23) 

and 

a  =  tan-1  (^2/^1) 

(24) 
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as  the  parameters  to  be  estimated.  The  log  likelihood  function  is  of  the  form 


L 

log m  =  -Xlog  |/„(«)|  -  £  YHAWl-Y,, 


*=1 


(25) 


where  the  form  taken  by  the  spectral  density  matrix  of  the  data  can  be  either  (8)  or  (12) 
depending  on  whether  we  have  the  perfect  or  general  correlation  models.  Suppose,  for 
the  moment,  that  we  only  want  estimators  for  the  velocity  and  azimuth.  Large  sample 
likelihood  theory  implies  that,  for  a  true  value  of  9,  the  distribution  of  9  is  approximately 
normal  with  mean  9  and  covariance  matrix 

,  f  ..a’logiwr1 

C00(»)  =  |--E  gggg,  }  (26) 

Having  done  this,  note  that  the  velocity  and  azimuth  axe  nonlinear  functions,  say  h(0)  = 
(c,<*y  =  (hi(0),h2(9))'  and  the  delta  method  implies  that  the  function  is  approximately 
normal  with  mean  h(#o)  and 


cov(b(t>)  =  (f£)  <*>”(«)  (f£) 
For  the  velocity  and  azimuth  functions,  note  that 

sh  =  j_(  \ 

90  ||s|3  Vii«ii^ 


(27) 


(28) 


The  following  sections  discuss  maximum  likelihood  estimation  and  derive  the  limiting 
distribution  of  the  maximum  likelihood  estimators  for  velocity  and  azimuth.  We  also  derive 
the  likelihood  ratio  detectors  for  the  signal  and  its  distribution  for  the  perfectly  correlated 
case. 

3.1  Maximum  Likelihood:  The  Perfectly  Correlated  Signal 

For  the  perfectly  correlated  case,  with  covariance  matrix  (11),  we  may  write  the  log 
likelihood  function  (25)  as 


logX(0)  =  —L  log 


P,  1  Av.v  1  EL  lg-(»)Y<P 

fTTT- jrLY‘Y‘  + 


(29) 


(N  +  r)  '  Pn  N  +  r 

This  is  seen  to  be  a  monotone  function  of  the  beam  power  which  will  be  proportional  to 

(30) 


*(*)  =  £!  g*(*)Yd2. 


<=1 


34 


Hence,  to  maximize  the  log  likelihood,  it  will  be  sufficient  to  maximize  the  beam  power 
(30)  over  9. 

Suppose  that  we  look  for  the  liklihood  ratio  criterion  for  testing  the  presence  or  absence 
of  the  signal  St-  Then,  for  St  =  0,  we  will  obtain  a  test  statistic  that  is  a  monotone  function 
of  the  beam  power  in  (30).  For  any  9 ,  under  the  hypothesis  P3  =  0,  the  distribution  of 
B(9)  is  proportional  to  a  chi-squared  distribution  with  2 L  degrees  of  freedom,  i.e. 


2  B{9) 
NPn 


~  X2(2 L). 


(31) 


Under  the  alternative  hypthesis,  assuming  the  wavenumber  vector  is  evaluated  at  the 
correct  9 ,  we  have 


2  B{9) 
NPn 


(32) 


If  9q  is  the  model  value  and  we  use  the  beam  at  9 ,  the  distribution  of  the  test  statistic  is 


2  B(9) 
NPn 


(l  +  d(9, 0o)p“ )} 

■*; n 


(33) 


where 


N  N 


d(9, 90)  =  N  1  cos(27r(ri  -  r*)'(^  ~  «o)] 

j=l  Jk=l 


(34) 


and  the  detection  probability  is  a  function  of  the  offset  between  9  and  90. 

The  uncertainty  of  the  maximum  likelihood  estimators  for  velocity  and  azimuth  are 
evaluated  by  using  (26)-(28)  in  conjunction  with  the  log  likelihood  (29)  and  we  note  that 
the  covariance  matrix  of  9  simplifies,  since 


a2  log  L(9) 
d9d9’ 


(2^)2 

Pn(N  +  r) 


L 

Yj  Y  exP{2 **(ri  -  rk),9}YjiYCt( Tj  -  rjk)(ry  -  rk)' . 

j,k  /=1 


Hence, 


—E 


l 


a2  log  j(gn  2  l  p, 

dddd'  J  [  }  (N  +  r)Pn 


Yw 


-  rir k  -  rkr'j  +  Tfcrjt) 


_  2(27r)2iV2L 
r(N  +  r) 


where 


R  = 


^  N 

N 

J=  1 


f)(r, 


(35) 


(36) 
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is  the  sample  covariance  matrix  of  the  array  of  the  array  coordinates.  It  follows  that  9  will 
be  approximately  normal  with  mean  zero  and  approximate  covariance  matrix 


cov(9) 


2(2tt)2  L  N 


K> 


where  r  is  the  inverse  of  the  signal-to-noise  ratio  (13).  Then,  defining  the  vectors  8  — 
(0l5  02y  and  9  =  (02>  — 0i)',  we  obtain  the  variance  estimators  for 


var  c 


1  r  /  r\  v2  8'R-l9 
2(27t)2  N  \  +  N )  ||0||6  L 


(38) 


and 


1  r  /  r  \  1  9'R~l9 

2t2rj*N  V+  NJ  ||0||4  L  ' 


(39) 


Evaluating  the  above  equations  at  9  will  produce  consistent  estimators  for  the  variances. 


3.1  Maximum  Likelihood:  The  General  Correlated  Signal 

For  the  signal  with  a  general  correlation  structure,  as  in  (1),  the  log  likelihood  has  the 
covariance  structure  given  by  (7)  and  we  write  (24)  in  the  form 

logL(0)  oc  — Llog  |/3 1  —  log|I/\r  +  -Pn/T1! 

-  i-  V  Y*eYt  Y  *tG{9)(lN  +  Pnf;1)  G*(8)Yt.  (40) 

Pn  (  Pn  l  '  ' 

Then,  noting  that  the  Hermitian  form  contains  the  matrix  C  in  (21),  we  may  write  the 
likelihood  ratio  detector  in  the  form 

i=i 

using  (18).  The  detector  does  not  seem  to  have  a  tractable  distribution  and  it  may  be 
better  to  compare  2logL(0)  -2  log  I,  with  a  xl  distribution  where  logL  denotes  the  value 
of  (40)  evalulated  at  fa  —  0. 

We  may  derive  variance  formulae  using  an  argument  similar  to  that  for  the  perfectly 
correlated  case.  First,  note  that 


^  log  L(9) 
8988' 


^{>|Y;WW)y,} 

-  *  w} 

^  j,k  t 


=  -  r»)(rJ  -  n)' 

"  i.t  < 
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Taking  expectations,  we  obtain 


= w2k  T,e*Mri  - ~  r*>' 

=  (2  *?±rD, 

i  n 

where  fh  denotes  the  jkth  element  of  the  signal  spectral  matrix  and 


D  =  5Z  “  rk)(ri  ~  r*)*' 


(42) 


(43) 


In  this  case  we  will  have  9  distributed  approximately  as  a  normal  random  variable  with 
mean  9q  and  covariance  matrix 


cov  9 


1  Pn 


D 


(2tt)2  L 

leading  to  estimated  variances  for  the  velocity  and  azimuth  of  the  form 


var  c 


1  „  i/2  9'D~l9 

Ai 


(2?r)2  "||*|| 


and 


var  a  ~ 


1  1  9'D~X9 

(27t)2  n\\9\\*  L 


(44) 


(45) 


(46) 
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APPENDIX:  Variance  of  Beamforming  Estimator  Under  Signed  Decorrelation 

For  the  original  estimated  wavenumber  in  the  decorrelated  signal  case,  we  have 

where 

D  =  r,-  -  ri)(r,-  -  r ,•)' 

with 

c  =  (In  +  Pnf:1)-' 

For  the  maximizer  of  the  beam  power,  say 

m = E  te*wY<i2- 

e=i 


say  i  we  have 


where 


and 


We  get 


and 


BB-AN^-^V-'WV-'y 
W=  E(*>~  ri)(r<'  -  */)'!!>■  lh 

+  (r«  “  rj)(ri'  -  *j)'fii'Pn  +  “  rj)(ri  ~  rj'YfPjPn 

+  S(r‘  “  r>)(r‘  “  ri)'Pn 


3,r,* 


F  =  X)(r«  “  ri)(r>  “  ri)7y 


«*J 


1  „  1  O'D^d 

var  a  »  —  ,0Pn 


var  as 


(2»)»-"||«||4  I 
1  1 

(27T)2  ||*||4  L 
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function  y  =  fsnewth(c,T,B,tw,sni,ax,deg,dc) 

%  function  y  =  snewth(c,T,B,tw,sni,xx) 

%  standard  error  in  degrees  of  optimum  beamed  estimate 
&  of  signal  with  imperfect  correlation 
%  fs  returns  matrix  of  coherences  between  sensors 
%  c:  velocity  km/sec,  T:  signal  period  seconds, 

%  k  equals  1/wavelength,  not  2*pi/wavelength 

%  ax  is  N  rows  by  2  columns,  x  and  y  coordinates  in  km 

%  B:  bandwidth,  tw:  time  window,  sni:  amplitude  S/N  on  each  element 

%  N:  number  of  sensors  in  array 

%  az  is  0  to  the  north  and  increases  clockwise, 

%  deg:  1/2  delta  azimuth  e.g.  5  degrees 
%  dc:  delta  velocity  in  km/sec,  e.g.  for  infrasound  0.01 5 

%  d  is  D/Ps  and  is  the  coherence-weighted  array  covariance  matrix  corresponding 
%  to  R;  Ps  is  the  signal  power  spectrum 


% 

%y=(180/pi)*(l/(2*pi))*((c*T)/sqrt(xx))*sqrt((l/(2*B*tw))*(l/(N*sniA2))*(l+(l/(N*sniA2))»; 


N=size(ax,l); 
fsm=fs(c,T,az,ax,deg,dc) ; 
cc=inv(eye(N)+(  1  /sni  A2)*(inv(fsm))) : 
d=zeros(2); 
for  id=l  :2; 
for  jd=l:2; 
for  j=l:N; 
for  k=l  :N; 

d(idjd)=d(idjd)+cc(j,k)*fsm(kj)*(ax(j,id)-ax(k,id))*(ax(jjd)-ax(kjd)); 

end 

end 

end 

end 


ko=l/(c*T);  kx=ko*sin(pi/2-az);  ky=ko*cos(pi/2-az); 
theta=[ky  -kx]’; 
kidk=theta’  *inv(d)*theta; 

y=(180/pi)*(l/(2*pi))*(l/sqrt(sniA2*B*tw))*(l/koA2)*sqrt(kidk); 
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function  y  =  bfsnewth(c,T,B,tw,sni,ax,az,deg,dc) 

%  standard  error  in  degrees  of  beamed  estimate 
&  of  signal  with  imperfect  correlation 
%  fs  returns  matrix  of  coherences  between  sensors 
%  c:  velocity  km/sec,  T:  signal  period  seconds, 

%  k  equals  1 /wavelength,  not  2*pi/wavelength 
%  ax  is  N  rows  by  2  columns,  x  and  y  coordinates  in  km 
%  B:  bandwidth,  tw:  time  window,  sni:  *amplitude*  S/N  on  each  element 
%  (in  text  snr  is  *power*  S/N,  therefore  sniA2  in  formulas  below) 

%  N:  number  of  sensors  in  array 
%  az  is  0  to  the  north  and  increases  clockwise, 

%  deg:  1/2  delta  azimuth  e.g.  5  degrees 
%  dc:  delta  velocity  in  km/sec,  e.g.  for  infrasound  0.01 5 
%  v  is  V/(Pn*sniA2)  is  the  coherence- weighted  array  covariance  matrix 
%  w  is  W/Pn**2,  note  in  final  formula  sniA4  instead  of  sniA2  in 
%  optimum  formula  "because"  of  sniA4  in  w 

%  Note  for  perfect  correlation: 

%  y=(180/pi)*(l/(2*pi))*sqrt((l/(2*B*tw))*(l/(N*sniA2)) 

%  *(l+(l/N*sniA2))))*((c*T)/sqrt(xx)); 


N=size(ax,l); 

fsm=fs(c,T,az,ax,deg,dc) ; 

v=zeros(2); 

w=zeros(2); 

for  id=l  :2; 

for  jd=l:2; 

for  i=l:N; 

for  j=l:N; 

v(id,jd)=v(idjd)+fsm(ij)*(ax(i,id)*(ax(j,id)-ax(jjd»; 

end 

end 

end 

end 

for  id=l  :2; 
for  jd=l:2; 
for  i=l:N; 
for  j=l:N; 

w(idjd)=w(idjd)+  (ax(i,id)-ax(j,id)*(ax(ijd)-ax(jjd)); 

end 

end 

end 

end 
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forid=l:2; 
for  jd=l:2; 
for  i=l:N; 
forj=l:N; 
for  jp=l:N; 

w(id,jd)=w(id,jd)+fsm(jp,j)*sniA2*(ax(i,id)-ax(j,id))*(ax(i,jd)-ax(jpjd)); 

end 

end 

end 

end 

end 

for  id=l  :2; 
for  jd=l:2; 
for  i=l:N; 
for  j=l:N; 
for  ip=l  :N; 

w(id,jd)=w(id,jd)+fsm(i,ip)*sniA2*(ax(i,id)-ax(j,id))*(ax(ip,jd)-ax(jjd)); 

end 

end 

end 

end 

end 

for  id=l  :2; 
for  jd=l:2; 
for  i=l:N; 
for  j=l:N; 
for  ip=l:N; 
forjp=l:N; 

w(idjd)=w(idjd)+fsm(i,ip)*sniA2*fsm(jpj)*sniA2*(ax(i,id)-ax(j,id))*(ax(ipjd)-ax(jp,jd)); 

end 

end 

end 

end 

end 

end 

ko=l/(c*T);  kx=ko*sin(pi/2-az);  ky=ko*cos(pi/2-az); 
theta=[ky  -kx]’; 

kidk=theta’  *  inv(v)  *  w*  inv(  v)  *  theta; 

y=(180/pi)*(l/(2*pi))*(l/sqrt(sniA4*B*tw))*(l/koA2)*sqrt(kidk); 
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function  fsm=fs(c,T,az,a,deg,dc) 


%  coherence  matrix  for  array  azimuth  estimation 
%  a  is  array  element  location  x,y  matrix 
%  a  is  N  rows  by  2  columns,  x  and  y  coordinates  in  km 

%  c:  velocity  km/sec,  T:  signal  period  seconds,  az:  event  azimuth 
%  az  is  0  to  the  north  and  increases  clockwise, 

%  xl,x2,yl,y2:  coordinates  of  sensors, 

%  deg:  1/2  delta  azimuth  e.g.  5  degrees 
%  dc:  delta  velocity  in  km/sec,  e.g.  for  infrasound  0 

N=size(a,l); 
for  i=l  :N 
forj=l:N 

xl  =a(i,  1  );x2=a(j,  1  );y  1  =a(i,2);  y2=a0',2); 
fsm(i,j)=(gams(c,T,az,x  1  ,x2,y  1  ,y2,deg,dc»; 
end 
end 


function  z  =  gams(c,T,az,xl,x2,yl,y2,deg,dc) 

%  function  z  =  gams(c,T,az,xl,x2,yl,y2,deg,dc) 

%  squared  coherence  of  signal  propagating  between  two  stations 
%  Mack  and  Flinn  Geophysical  Journal  (1971),  26,  p  263,  eqn  (3) 

%  c:  velocity  km/sec,  T:  signal  period  seconds,  az:  event  azimuth  in  radians, 
%  az  is  0  to  the  north  and  increases  clockwise 
%  az  in  radians,  e.g.  30  degrees  is  pi/6  (modified  so  az  in  degrees) 

%  thr  is  0  to  east  and  increases  counter-clockwise 

%  xl,x2,yl,y2:  coordinates  of  sensors,  deg:  1/2  delta  azimuth  e.g.  5  degrees 
%  dc:  delta  velocity  in  km/sec,  e.g.  for  infrasound  0.015 

az=(az/90)*pi/2; 
ko=l/(c*T); 
dth  =  pi*deg/180; 

%  dth  is  1/2  delta  azimuth  in  radians 

dk  =  l/(c*T)  -l/((c+dc)*T); 

thr=atan2(y2-yl,x2-xl); 
thd=(pi/2  -  az)  -  thr; 

%  thd  is  angle  from  r  to  k,  positive  counter-clockwise 

r=sqrt((y2-yl  )A2+(x2-xl  )A2); 
y=r*cos(thd);x=r*sin(thd); 


if  x  ==  0;  x=eps;  end; 
if  y  =  0;  y=eps;  end; 


argl=2*pi*ko*x*sin(dth);  arg2=2*pi*dk*y; 


zr=(sin(argl  )/argl  )*(sin(arg2)/ arg2); 
z=zrA2 
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